import pandas as pd
import numpy as np
import geopandas as gpd
import seaborn as sns
import matplotlib.pyplot as plt
from itables import init_notebook_mode
import itables
import warnings
warnings.filterwarnings("ignore", "use_inf_as_na")
init_notebook_mode(all_interactive=True)
import plotly.express as px
df = pd.read_excel('data/who_aap_2021_v9_11august2022.xlsx', sheet_name='AAP_2022_city_v9')
df.columns
Index(['WHO Region', 'ISO3', 'WHO Country Name', 'City or Locality',
'Measurement Year', 'PM2.5 (μg/m3)', 'PM10 (μg/m3)', 'NO2 (μg/m3)',
'PM25 temporal coverage (%)', 'PM10 temporal coverage (%)',
'NO2 temporal coverage (%)', 'Reference',
'Number and type of monitoring stations', 'Version of the database',
'Status'],
dtype='object')
itables.show(df)
| WHO Region | ISO3 | WHO Country Name | City or Locality | Measurement Year | PM2.5 (μg/m3) | PM10 (μg/m3) | NO2 (μg/m3) | PM25 temporal coverage (%) | PM10 temporal coverage (%) | NO2 temporal coverage (%) | Reference | Number and type of monitoring stations | Version of the database | Status |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Loading... (need help?) |
Biblioteka ITables 2.0 stanowi użyteczne narzędzie - można przeglądać dane łatwiej niż za pomocą scrollowania i zadecydować ile chcemy ich widzieć na raz. Nie jest ona bezwzględnie potrzebna do analizy danych, ale na pewno stanowi ułatwienie pod względem czytelności i przejrzystości.
df.info()
<class 'pandas.core.frame.DataFrame'> RangeIndex: 32191 entries, 0 to 32190 Data columns (total 15 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 WHO Region 32190 non-null object 1 ISO3 32191 non-null object 2 WHO Country Name 32191 non-null object 3 City or Locality 32191 non-null object 4 Measurement Year 32191 non-null int64 5 PM2.5 (μg/m3) 15048 non-null float64 6 PM10 (μg/m3) 21109 non-null float64 7 NO2 (μg/m3) 22200 non-null float64 8 PM25 temporal coverage (%) 7275 non-null float64 9 PM10 temporal coverage (%) 5381 non-null float64 10 NO2 temporal coverage (%) 19890 non-null float64 11 Reference 32186 non-null object 12 Number and type of monitoring stations 8758 non-null object 13 Version of the database 32191 non-null int64 14 Status 0 non-null float64 dtypes: float64(7), int64(2), object(6) memory usage: 3.7+ MB
df.describe()
| Measurement Year | PM2.5 (μg/m3) | PM10 (μg/m3) | NO2 (μg/m3) | PM25 temporal coverage (%) | PM10 temporal coverage (%) | NO2 temporal coverage (%) | Version of the database | Status | |
|---|---|---|---|---|---|---|---|---|---|
| Loading... (need help?) |
Mimo widocznych wartości odstających zdecydowano o nieusuwaniu ich, ponieważ mogłoby to uniemożliwić zaobserwowanie ważnych wniosków - zanieczyszczenia mogą być na pewnych terenach ekstremalne.
print("Ilość wszystkich wierszy:", len(df))
df.isna().sum()
Ilość wszystkich wierszy: 32191
| 0 | |
|---|---|
| Loading... (need help?) |
Kolumna Status składa się tylko z NaN - można ją bezpiecznie usunąć. Kolumny WHO Region i Reference mają bardzo małą ilość NaN (<5%), te wiersze można wypełnić odpowiednimi informacjami (WHO Region - European Region dla Liechtensteinu) lub usunąć (Reference - brak informacji o tym, jak je wypełnić). Inne kolumny mają wiele wartości NaN (około > 10 000), więc zostaną pozostawione bez zmian. Są to najważniejsze kolumny dotyczące parametrów PM2.5, PM10 oraz NO2, więc można domyślać się, że istnieją braki w odczytywanych wartościach na czujnikach. Niestety sporo danych dotyczących pokrycia (części roku, w której zbierane są dane) również brakuje.
df.drop(columns=['Status'], inplace=True)
print(df["WHO Region"].unique())
df.loc[df["WHO Region"].isna(), "WHO Region"] = 'European Region'
['Eastern Mediterranean Region' 'European Region' 'Region of the Americas' 'Western Pacific Region' 'South East Asia Region' 'African Region' nan]
df.loc[df["Reference"].isna(), :] # ISO3 = QAT
| WHO Region | ISO3 | WHO Country Name | City or Locality | Measurement Year | PM2.5 (μg/m3) | PM10 (μg/m3) | NO2 (μg/m3) | PM25 temporal coverage (%) | PM10 temporal coverage (%) | NO2 temporal coverage (%) | Reference | Number and type of monitoring stations | Version of the database | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Loading... (need help?) |
df.loc[df["ISO3"]=="QAT", :]
| WHO Region | ISO3 | WHO Country Name | City or Locality | Measurement Year | PM2.5 (μg/m3) | PM10 (μg/m3) | NO2 (μg/m3) | PM25 temporal coverage (%) | PM10 temporal coverage (%) | NO2 temporal coverage (%) | Reference | Number and type of monitoring stations | Version of the database | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Loading... (need help?) |
df.dropna(subset=["Reference"],inplace=True)
# Double-check
print("Number of rows in the dataframe:", len(df))
df.isna().sum()
Number of rows in the dataframe: 32186
| 0 | |
|---|---|
| Loading... (need help?) |
Dodatkowa selekcja cech - zdecydowano odrzucić kolumnę Version Of The Database, ponieważ nie wnosi ona nic merytorycznego do analizy.
df.drop(columns=["Version of the database"], inplace=True)
Poniżej przedstawiono podstawowe wizualizacje danych (pairplot, histplot, heatmap), aby lepiej zrozumieć ich charakter przed dalszą analizą. Przedstawiono również jak wyglądają zanieczyszczenia dla poszczególnych lat i jak wyglądał coverage.
sns.pairplot(df)
plt.show()
Na pairplotcie widać relacje między zmiennymi oraz histogram każdej z nich. Zależność liniową widać pomiędzy PM2.5 oraz PM10. NO2 również posiada w jakimś stopniu zależność do PM2.5 jak i PM10, ale tego typu zależnośći zostaną zbadane dokładniej w dalszej kolejnośći. Histogramom przyjrzyjmy się dokładniej poniżej.
numeric_columns = df.select_dtypes(include=['float64', 'int64'])
cols = numeric_columns.columns.to_list()
fig, ax = plt.subplots(2,4,figsize=(20,10))
ax = ax.flatten()
for i in range(len(cols)):
sns.histplot(df[cols[i]], ax=ax[i])
ax[i].set(title=f'Histplot {cols[i]}', ylabel='counts')
plt.show()
Pierwszym wnioskiem jest obserwacja na temat czasu kiedy były zbierane dane - sensowna ilość obserwacji zaczyna pojawiać się dopiero od 2010 roku (być może wtedy nastąpił skok technologiczny dotyczący czujników lub pojawiła się większa ich dostępność na świecie). Obserwacje sprzed 2010 zostaną usunięte. Warto zauważyć, że w latach 2011-12 zanotowano spadek obserwacji a następnie nagły skok aż do 2020 kiedy wybuchła pandemia (spadek obserwacji może wynikać właśnie z niej i powodanego przez nią kryzysu na świecie - były wtedy ważniejsze sprawy do monitorowania niż jakość powietrza). Jeśli chodzi o zmiennie dotyczące PM2.5, PM10 oraz NO2 mają one rozkłady prawoskośne, natomiast wykresy dotyczące pokrycia mają rozkłady lewoskośne. To dobrze, ponieważ im mniejsze zanieczyszczenie tym lepiej dla świata, a im większe pokrycie tym lepiej, bo dane są dokładniejsze i pozwalają na bardziej wiarygodne analizy.
df2010 = df[df['Measurement Year']>=2010].reset_index(drop=True)
print(len(df), len(df2010)) # niewielka róznica około 20 rekordów
32186 32165
numeric_columns = df2010.select_dtypes(include=['float64', 'int64'])
correlation = numeric_columns.corr()
plt.figure(figsize=(10, 8))
sns.heatmap(correlation, cmap='coolwarm', fmt=".2f", vmin=-1, vmax=1)
plt.title('Heatmapa korelacji zmiennych numerycznych')
plt.show()
Z heatmapy wypływa podobny wniosek jak z pairplotu - występuje związek między PM2.5 oraz PM10, a w mniejszym stopniu między NO2 oraz PM2.5 jak i NO2 oraz PM10. Podobnie dla tych samych par wygląda sytuacja dla zmiennych dotyczących pokrycia.
Poniżej przedstawiono jak wygląda zanieczyszczenie na przestrzeni lat.
columns = ['PM2.5 (μg/m3)', 'PM10 (μg/m3)', 'NO2 (μg/m3)']
plt.figure(figsize=(10, 6))
for column in columns:
sns.lineplot(data=df2010, x='Measurement Year', y=column, label=column, errorbar=None)
plt.title('Szeregi czasowe dla PM2.5, PM10 i NO2 dla całych danych')
plt.xlabel('Rok')
plt.ylabel('Wartość')
plt.legend(title='Parametr')
plt.tight_layout()
plt.show()
Największe zanieczyszczenie powoduje PM10, natomiast PM2.5 oraz NO2 wypadają dosyć podobnie. Można zauważyć, że nie ma wielu wspólnych tendencji dla wszystkich 3 parametrów na raz - ewentualnie stosunkowa stabilność w latach 2014-2017.
columns_cov = ['PM25 temporal coverage (%)', 'PM10 temporal coverage (%)', 'NO2 temporal coverage (%)']
plt.figure(figsize=(10, 6))
for column in columns_cov:
sns.lineplot(data=df2010, x='Measurement Year', y=column, label=column, errorbar=None)
plt.title('Szeregi czasowe dla pokrycia czasowego PM2.5, PM10 i NO2 dla całych danych')
plt.xlabel('Rok')
plt.ylabel('Pokrycie czasowe (%)')
plt.legend(title='Parametr')
plt.tight_layout()
plt.show()
Dla NO2 pokrycie wygląda najstabilniej - dane po podejrzanych latach 2011-12 utrzymywał się aż do pandemii na stałym wysokim poziomie, dane były zbierane przez więcej niż 90% czasu w roku. Najmniej stabilnie przedstawia się PM2.5, może to wynikać z jakości czujników lub problemów technicznych.
len(df2010['WHO Region'].unique()) # 6 regionów - wizualizacje będą czytelne
6
x_ticks = [x.replace(" ", "\n") for x in df2010["WHO Region"].unique()]
sns.histplot(data=df2010, x="WHO Region",bins=20, alpha=0.5)
plt.xticks(ticks=np.arange(len(x_ticks)),labels=x_ticks)
plt.title(f'Distribution of WHO Region values')
plt.show()
W regionie europejskim zebrano zdecydowanie najwięcej danych, należy pamiętać o tym podczas dalszych analiz - gdy patrzymy na dane całościowo może wystąpić "europejski bias". Poniżej przedstawionio wykresy pudełkowe zmiennych PM2.5, PM10, NO2 z podziałem na regiony oraz rozkładem danej zmiennej dla odniesienia.
fig, ax = plt.subplots(3,2, figsize=(10,15))
for i, column in enumerate(columns):
sns.boxplot(data=df2010, y=column, x='WHO Region', ax=ax[i,0])
ax[i,0].set_xticks(ticks=np.arange(len(x_ticks)), labels=x_ticks)
ax[i,0].set_title(f'Rozkład {column} z podziałem na regiony')
sns.histplot(data=df2010, x=column, kde=True, bins=20, ax=ax[i,1], alpha=0.5)
ax[i,1].set_title(f'Rozkład {column}')
plt.tight_layout()
plt.show()
Z wizualizacji można zauważyć, że większość obserwacji skupia się wokół mniejszych wartości, ale widać sporo outlierów w górnej części wykresów pudełkowych.
fig, ax = plt.subplots(3, 1, figsize=(15, 15))
for i, column in enumerate(columns):
sns.lineplot(data=df2010, x='Measurement Year', y=column, hue='WHO Region', ax=ax[i], errorbar=None, palette='tab20')
ax[i].set_title(f'Szeregi czasowe dla {column} z podziałem na regiony')
ax[i].set_ylabel(column)
ax[i].set_xlabel('Rok')
ax[i].legend(bbox_to_anchor=(1.05, 1), loc='upper left')
plt.tight_layout()
plt.show()
Napawającym optymizmem wnioskiem jest to, że z czasem zanieczyszczenie ma tendencje do spadku (wszystkie 3 parametry).
len(df2010['WHO Country Name'].unique())
118
Krajów jest za dużo, aby czytelnie zwizualizować ich cechy. Postanowiono wybrać po 2 kraje dla każdego regionu z największą ilością rekordów, aby były reprezentatywne.
def top_countries_with_counts(group):
counts = group['WHO Country Name'].value_counts().nlargest(2)
return list(zip(counts.index, counts.values))
top_countries_by_region = df2010.groupby('WHO Region').apply(top_countries_with_counts).reset_index(name='Top Countries')
top_countries = [country[0] for sublist in top_countries_by_region['Top Countries'] for country in sublist]
top_countries = list(set(top_countries))
df_countries = df2010[df2010['WHO Country Name'].isin(top_countries)]
print(len(df_countries))
top_countries_by_region
15492
| WHO Region | Top Countries |
|---|---|
| Loading... (need help?) |
Warto zauważyć, że dla Afryki mamy bardzo mało obserwacji w porównaniu do innych regionów - aby o tym pamiętać pozostawiono wybrane kraje do dalszej analizy. Poniżej rozkład wartości dla wybranych państw.
x_ticks = [x.replace(" ", "\n") for x in df_countries["WHO Country Name"].unique()]
fig, ax = plt.subplots(1,1,figsize=(12,5))
sns.histplot(ax=ax,data=df_countries, x="WHO Country Name",bins=20, alpha=0.5)
plt.xticks(ticks=np.arange(len(x_ticks)),labels=x_ticks, rotation=45)
plt.title(f'Distribution of WHO Country Name values')
plt.show()
columns = ['PM2.5 (μg/m3)', 'PM10 (μg/m3)', 'NO2 (μg/m3)']
fig, ax = plt.subplots(3,1, figsize=(10,15))
for i, column in enumerate(columns):
sns.boxplot(data=df_countries, y=column, x='WHO Country Name', ax=ax[i])
ax[i].set_xticks(ticks=np.arange(len(x_ticks)), labels=x_ticks)
ax[i].set_title(f'Rozkład {column} z podziałem na regiony')
ax[i].set_xticklabels(x_ticks, rotation = 50)
plt.tight_layout()
plt.show()
Podobnie jak przy analizie regionów widać, że większość wartości skupia się w dolnej części boxplotów. Natomiast na tym poziomie granularności widać rozbieżności w ilości pomiarów dla poszczególnych regionów.
fig, ax = plt.subplots(3, 1, figsize=(15, 15))
for i, column in enumerate(columns):
sns.lineplot(data=df_countries, x='Measurement Year', y=column, hue='WHO Country Name', ax=ax[i], errorbar=None, palette='tab20')
ax[i].set_title(f'Szeregi czasowe dla {column} z podziałem na reprezentatywne kraje')
ax[i].set_ylabel(column)
ax[i].set_xlabel('Rok')
ax[i].legend(bbox_to_anchor=(1.05, 1), loc='upper left')
plt.tight_layout()
plt.show()
len(df2010['City or Locality'].unique())
6871
Analogicznie jak dla miast zostaną wybrane tylko reprezentatywne państwa w oparciu o regiony.
def top_cities_with_counts(group):
counts = group.groupby(['City or Locality', 'WHO Country Name']).size().nlargest(2)
return list(zip(counts.index.get_level_values(0), counts.index.get_level_values(1), counts.values))
top_cities_by_region = df2010.groupby('WHO Region').apply(top_cities_with_counts).reset_index(name='Top Cities')
top_cities = [city[0] for sublist in top_cities_by_region['Top Cities'] for city in sublist]
top_cities = list(set(top_cities))
df_cities = df2010[df2010['City or Locality'].isin(top_cities)]
print(len(df_cities))
top_cities_by_region
146
| WHO Region | Top Cities |
|---|---|
| Loading... (need help?) |
Jako, że dane są zbierane co rok, a widać wartości np. 17 to widać, że coś jest nie w porządku, postanowiono przyjrzeć się danym.
df_cities
| WHO Region | ISO3 | WHO Country Name | City or Locality | Measurement Year | PM2.5 (μg/m3) | PM10 (μg/m3) | NO2 (μg/m3) | PM25 temporal coverage (%) | PM10 temporal coverage (%) | NO2 temporal coverage (%) | Reference | Number and type of monitoring stations | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Loading... (need help?) |
Dane zawierają duplikaty z różnicą jedynie w kolumnie Number and type of monitoring stations. Jako, że kluczowe parametry są te same postanowiono pozbyć się tej kolumny a następnie usunąć duplikaty.
df_unique = df2010.drop(columns=['Number and type of monitoring stations'])
df_unique = df_unique.drop_duplicates()
print(len(df2010),len(df_unique))
32165 32069
Widać, że zniknęło około 100 rekordów - nie ma to większego wpływu na poprzednie analizy, więc nie wprowadzano żadnych korekt.
top_cities_by_region = df_unique.groupby('WHO Region').apply(top_cities_with_counts).reset_index(name='Top Cities')
top_cities = [city[0] for sublist in top_cities_by_region['Top Cities'] for city in sublist]
top_cities = list(set(top_cities))
df_cities = df_unique[df_unique['City or Locality'].isin(top_cities)]
print(len(df_cities))
top_cities_by_region
130
| WHO Region | Top Cities |
|---|---|
| Loading... (need help?) |
fig, ax = plt.subplots(3, 1, figsize=(15, 15))
for i, column in enumerate(columns):
sns.lineplot(data=df_cities, x='Measurement Year', y=column, hue='City or Locality', ax=ax[i], errorbar=None, palette='tab20')
handles, labels = ax[i].get_legend_handles_labels()
new_labels = [f"{label} ({df_cities[df_cities['City or Locality'] == label]['WHO Country Name'].iloc[0]})" for label in labels]
ax[i].legend(handles, new_labels, bbox_to_anchor=(1.05, 1), loc='upper left')
ax[i].set_title(f'Szeregi czasowe dla {column} z podziałem na reprezentatywne miasta')
ax[i].set_ylabel(column)
ax[i].set_xlabel('Rok')
plt.tight_layout()
plt.show()
mean_pollution = df_unique.groupby(['WHO Country Name','ISO3'])[['PM2.5 (μg/m3)', 'PM10 (μg/m3)', 'NO2 (μg/m3)']].mean().reset_index()
for col in columns:
fig = px.choropleth(
mean_pollution,
locations='ISO3',
locationmode='ISO-3',
color=col,
hover_name='WHO Country Name',
projection='natural earth',
title=f'Średnie wartości {col} dla całości danych (lata 2010-2021)'
)
fig.show()
Na pytania wylistowane na UPELu częsciowo odpowiedzinao w komentarzach i analizach powyżej. Poniżej jednak znajduje się podsumowanie i wypunktowane najważniejsze wnioski:
Związki przyczynowo-skutkowe:
Kontekst ML: